################################################################################################################
#'@title  Data-Driven Nonparametric Shape Restriction and Parametric Model Specification Testing using Binscatter
#'@description \code{binstest} implements binscatter-based hypothesis testing procedures for parametric functional
#'             forms of and nonparametric shape restrictions on the regression function of interest, following the results
#'             in \href{https://arxiv.org/abs/1902.09608}{Cattaneo, Crump, Farrell and Feng (2021a)}.
#'             If the binning scheme is not set by the user,
#'             the companion function \code{\link{binsregselect}} is used to implement binscatter in a
#'             data-driven way and inference procedures are based on robust bias correction.
#'             Binned scatter plots based on different methods can be constructed using the companion functions \code{\link{binsreg}},
#'             \code{\link{binsqreg}} or \code{\link{binsglm}}.
#'@param  y outcome variable. A vector.
#'@param  x independent variable of interest. A vector.
#'@param  w control variables. A matrix, a vector or a \code{\link{formula}}.
#'@param  data an optional data frame containing variables used in the model.
#'@param  estmethod estimation method. The default is \code{estmethod="reg"} for tests based on binscatter least squares regression. Other options are \code{"qreg"} for quantile regression and \code{"glm"} for generalized linear regression. If \code{estmethod="glm"}, the option \code{family} must be specified.
#'@param  family a description of the error distribution and link function to be used in the generalized linear model when \code{estmethod="glm"}. (See \code{\link{family}} for details of family functions.)
#'@param  quantile the quantile to be estimated. A number strictly between 0 and 1.
#'@param  deriv  derivative order of the regression function for estimation, testing and plotting.
#'               The default is \code{deriv=0}, which corresponds to the function itself.
#'@param  at value of \code{w} at which the estimated function is evaluated.  The default is \code{at="mean"}, which corresponds to
#'           the mean of \code{w}. Other options are: \code{at="median"} for the median of \code{w}, \code{at="zero"} for a vector of zeros.
#'           \code{at} can also be a vector of the same length as the number of columns of \code{w} (if \code{w} is a matrix) or a data frame containing the same variables as specified in \code{w} (when
#'           \code{data} is specified). Note that when \code{at="mean"} or \code{at="median"}, all factor variables (if specified) are excluded from the evaluation (set as zero).
#'@param  nolink if true, the function within the inverse link function is reported instead of the conditional mean function for the outcome.
#'@param  testmodel a vector. \code{testmodel=c(p,s)} sets a piecewise polynomial of degree \code{p} with \code{s}
#'                  smoothness constraints for parametric model specification testing.  The default is
#'                  \code{testmodel=c(3,3)}, which corresponds to a cubic B-spline estimate of the regression
#'                  function of interest for testing against the fitting from a parametric model specification.
#'@param  testmodelparfit a data frame or matrix which contains the evaluation grid and fitted values of the model(s) to be
#'                        tested against.  The column contains a series of evaluation points
#'                        at which the binscatter model and the parametric model of interest are compared with
#'                        each other.  Each parametric model is represented by other columns, which must
#'                        contain the fitted values at the corresponding evaluation points.
#'@param  testmodelpoly degree of a global polynomial model to be tested against.
#'@param  testshape a vector. \code{testshape=c(p,s)} sets a piecewise polynomial of degree \code{p} with \code{s}
#'                  smoothness constraints for nonparametric shape restriction testing. The default is
#'                  \code{testshape=c(3,3)}, which corresponds to a cubic B-spline estimate of the regression
#'                  function of interest for one-sided or two-sided testing.
#'@param  testshapel a vector of null boundary values for hypothesis testing. Each number \code{a} in the vector
#'                   corresponds to one boundary of a one-sided hypothesis test to the left of the form
#'                   \code{H0: sup_x mu(x)<=a}.
#'@param  testshaper a vector of null boundary values for hypothesis testing. Each number \code{a} in the vector
#'                   corresponds to one boundary of a one-sided hypothesis test to the right of the form
#'                   \code{H0: inf_x mu(x)>=a}.
#'@param  testshape2 a vector of null boundary values for hypothesis testing. Each number \code{a} in the vector
#'                   corresponds to one boundary of a two-sided hypothesis test ofthe form
#'                   \code{H0: sup_x |mu(x)-a|=0}.
#'@param  lp an Lp metric used for (two-sided) parametric model specification testing and/or shape restriction testing. The default is \code{lp=Inf}, which
#'           corresponds to the sup-norm of the t-statistic. Other options are \code{lp=q} for a positive integer \code{q}.
#'@param  bins A vector. Degree and smoothness for bin selection. The default is \code{bins=c(2,2)}, which corresponds to a quadratic spline estimate.
#'@param  nbins number of bins for partitioning/binning of \code{x}.  If not specified, the number of bins is
#'              selected via the companion function \code{binsregselect} in a data-driven, optimal way whenever possible.
#'@param  binspos position of binning knots. The default is \code{binspos="qs"}, which corresponds to quantile-spaced
#'                binning (canonical binscatter).  The other options are \code{"es"} for evenly-spaced binning, or
#'                a vector for manual specification of the positions of inner knots (which must be within the range of
#'                \code{x}).
#'@param  binsmethod method for data-driven selection of the number of bins. The default is \code{binsmethod="dpi"},
#'                   which corresponds to the IMSE-optimal direct plug-in rule.  The other option is: \code{"rot"}
#'                   for rule of thumb implementation.
#'@param  nbinsrot initial number of bins value used to construct the DPI number of bins selector.
#'                 If not specified, the data-driven ROT selector is used instead.
#'@param  randcut upper bound on a uniformly distributed variable used to draw a subsample for bins selection.
#'                Observations for which \code{runif()<=#} are used. # must be between 0 and 1.
#'@param  nsims number of random draws for hypothesis testing. The default is
#'              \code{nsims=500}, which corresponds to 500 draws from a standard Gaussian random vector of size
#'              \code{[(p+1)*J - (J-1)*s]}.
#'@param  simsgrid number of evaluation points of an evenly-spaced grid within each bin used for evaluation of
#'                 the supremum (infimum or Lp metric) operation needed to construct hypothesis testing
#'                 procedures. The default is \code{simsgrid=20}, which corresponds to 20 evenly-spaced
#'                 evaluation points within each bin for approximating the supremum (infimum or Lp metric) operator.
#'@param  simsseed seed for simulation.
#'@param  vce procedure to compute the variance-covariance matrix estimator. For least squares regression and generalized linear regression, the allowed options are the same as that for \code{\link{binsreg}} or \code{\link{binsqreg}}.
#'            For quantile regression, the allowed options are the same as that for \code{\link{binsqreg}}.
#'@param  cluster cluster ID. Used for compute cluster-robust standard errors.
#'@param  asyvar  If true, the standard error of the nonparametric component is computed and the uncertainty related to control
#'                variables is omitted. Default is \code{asyvar=FALSE}, that is, the uncertainty related to control variables is taken into account.
#'@param  dfcheck adjustments for minimum effective sample size checks, which take into account number of unique
#'                values of \code{x} (i.e., number of mass points), number of clusters, and degrees of freedom of
#'                the different stat models considered. The default is \code{dfcheck=c(20, 30)}.
#'                See \href{https://nppackages.github.io/references/Cattaneo-Crump-Farrell-Feng_2021_Stata.pdf}{Cattaneo, Crump, Farrell and Feng (2021b)} for more details.
#'@param  masspoints how mass points in \code{x} are handled. Available options:
#'                   \itemize{
#'                   \item \code{"on"} all mass point and degrees of freedom checks are implemented. Default.
#'                   \item \code{"noadjust"} mass point checks and the corresponding effective sample size adjustments are omitted.
#'                   \item \code{"nolocalcheck"} within-bin mass point and degrees of freedom checks are omitted.
#'                   \item \code{"off"} "noadjust" and "nolocalcheck" are set simultaneously.
#'                   \item \code{"veryfew"} forces the function to proceed as if \code{x} has only a few number of mass points (i.e., distinct values).
#'                                          In other words, forces the function to proceed as if the mass point and degrees of freedom checks were failed.
#'                   }
#'@param  weights an optional vector of weights to be used in the fitting process. Should be \code{NULL} or
#'                a numeric vector. For more details, see \code{\link{lm}}.
#'@param  subset optional rule specifying a subset of observations to be used.
#'@param  numdist  Number of distinct for selection. Used to speed up computation.
#'@param  numclust Number of clusters for selection. Used to speed up computation.
#'@param  ...      optional arguments to control bootstrapping if \code{estmethod="qreg"} and \code{vce="boot"}. See \code{\link{boot.rq}}.
#'@return \item{\code{testshapeL}}{Results for \code{testshapel}, including: \code{testvalL}, null boundary values;
#'                                 \code{stat.shapeL}, test statistics; and \code{pval.shapeL}, p-value.}
#'        \item{\code{testshapeR}}{Results for \code{testshaper}, including: \code{testvalR}, null boundary values;
#'                                 \code{stat.shapeR}, test statistics; and \code{pval.shapeR}, p-value.}
#'        \item{\code{testshape2}}{Results for \code{testshape2}, including: \code{testval2}, null boundary values;
#'                                 \code{stat.shape2}, test statistics; and \code{pval.shape2}, p-value.}
#'        \item{\code{testpoly}}{Results for \code{testmodelpoly}, including: \code{testpoly}, the degree of global polynomial;
#'                               \code{stat.poly}, test statistic; \code{pval.poly}, p-value.}
#'        \item{\code{testmodel}}{Results for \code{testmodelparfit}, including: \code{stat.model}, test statistics;
#'                                \code{pval.model}, p-values.}
#'        \item{\code{opt}}{ A list containing options passed to the function, as well as total sample size \code{n},
#'                           number of distinct values \code{Ndist} in \code{x}, number of clusters \code{Nclust}, and
#'                           number of bins \code{nbins}.}
#'
#'@author
#' Matias D. Cattaneo, Princeton University, Princeton, NJ. \email{cattaneo@princeton.edu}.
#'
#' Richard K. Crump, Federal Reserve Bank of New York, New York, NY. \email{richard.crump@ny.frb.org}.
#'
#' Max H. Farrell, University of Chicago, Chicago, IL. \email{max.farrell@chicagobooth.edu}.
#'
#' Yingjie Feng (maintainer), Tsinghua University, Beijing, China. \email{fengyingjiepku@gmail.com}.
#'
#'@references
#' Cattaneo, M. D., R. K. Crump, M. H. Farrell, and Y. Feng. 2021a: \href{https://arxiv.org/abs/1902.09608}{On Binscatter}. Working Paper.
#'
#' Cattaneo, M. D., R. K. Crump, M. H. Farrell, and Y. Feng. 2021b: \href{https://arxiv.org/abs/1902.09615}{Binscatter Regressions}. Working Paper.
#'
#'@seealso \code{\link{binsreg}}, \code{\link{binsqreg}}, \code{\link{binsglm}}, \code{\link{binsregselect}}.
#'
#'@examples
#'  x <- runif(500); y <- sin(x)+rnorm(500)
#'  est <- binstest(y,x, testmodelpoly=1)
#'  summary(est)
#'@export

binstest <- function(y, x, w=NULL, data=NULL, estmethod="reg", family=gaussian(),
                     quantile=NULL, deriv=0, at=NULL, nolink=F,
                     testmodel=c(3,3), testmodelparfit=NULL, testmodelpoly=NULL,
                     testshape=c(3,3), testshapel=NULL,
                     testshaper=NULL, testshape2=NULL, lp=Inf,
                     bins=c(2,2), nbins=NULL, binspos="qs",
                     binsmethod="dpi", nbinsrot=NULL, randcut=NULL,
                     nsims=500, simsgrid=20, simsseed=NULL,
                     vce=NULL, cluster=NULL, asyvar=F,
                     dfcheck=c(20,30), masspoints="on", weights=NULL, subset=NULL,
                     numdist=NULL, numclust=NULL, ...) {

  # param for internal use
  qrot <- 2

  ####################
  ### prepare data ###
  ####################
  xname <- deparse(substitute(x))
  yname <- deparse(substitute(y))
  weightsname <- deparse(substitute(weights))
  subsetname  <- deparse(substitute(subset))
  clustername <- deparse(substitute(cluster))

  # extract y, x, w, weights, subset, if needed (w as a matrix, others as vectors)
  # generate design matrix for covariates W
  w.factor <- NULL
  if (is.null(data)) {
    if (is.data.frame(y)) y <- y[,1]
    if (is.data.frame(x)) x <- x[,1]
    if (!is.null(w)) {
      if (is.vector(w)) {
        w <- as.matrix(w)
        if (is.character(w)) {
          w <- model.matrix(~w)[,-1,drop=F]
          w.factor <- rep(T, ncol(w))
        }
      } else if (is.formula(w)) {
        w.model <- binsreg.model.mat(w)
        w <- w.model$design
        w.factor <- w.model$factor.colnum
      }
    }
  } else {
    if (xname %in% names(data)) x <- data[,xname]
    if (yname %in% names(data)) y <- data[,yname]
    if (weightsname != "NULL") if (weightsname %in% names(data)) {
      weights <- data[,weightsname]
    }
    if (subsetname != "NULL") if (subsetname %in% names(data))  {
      subset  <- data[,subsetname]
    }
    if (clustername != "NULL") if (clustername %in% names(data)) {
      cluster <- data[,clustername]
    }
    if (deparse(substitute(w))!="NULL") {
      if (is.formula(w)) {
        if (is.data.frame(at)) {
          if (ncol(data)!=ncol(at)) {
            data[setdiff(names(at), names(data))] <- NA
            at[setdiff(names(data), names(at))] <- NA
          }
          data <- rbind(data, at)
        }
        w.model <- binsreg.model.mat(w, data=data)
        w <- w.model$design
        w.factor <- w.model$factor.colnum
        if (is.data.frame(at)) {
          eval.w <- w[nrow(w),]
          w <- w[-nrow(w),, drop=F]
        }
      } else {
        if (deparse(substitute(w)) %in% names(data)) w <- data[,deparse(substitute(w))]
        w <- as.matrix(w)
        if (is.character(w)) {
          w <- model.matrix(~w)[,-1,drop=F]
          w.factor <- rep(T, ncol(w))
        }
      }
    }
  }

  if (!is.null(w)) nwvar <- ncol(w)
  else             nwvar <- 0

  # extract subset
  if (!is.null(subset)) {
    y <- y[subset]
    x <- x[subset]
    w <- w[subset, , drop = F]
    if (!is.null(cluster)) cluster <- cluster[subset]
    if (!is.null(weights)) weights <- weights[subset]
  }

  na.ok <- complete.cases(x) & complete.cases(y)
  if (!is.null(w)) na.ok <- na.ok & complete.cases(w)
  if (!is.null(weights)) na.ok <- na.ok & complete.cases(weights)
  if (!is.null(cluster)) na.ok <- na.ok & complete.cases(cluster)

  y <- y[na.ok]
  x <- x[na.ok]
  w <- w[na.ok, , drop = F]
  weights <- weights[na.ok]
  cluster <- cluster[na.ok]

  xmin <- min(x); xmax <- max(x)

  # evaluation point of w
  if (is.null(at))  at <- "mean"

  # change method name if needed
  if (!is.null(family)) if (family$family!="gaussian" | family$link!="identity") {
    estmethod <- "glm"
  }

  # extract family obj
  if (estmethod=="glm") {
    familyname <- family$family
    linkinv    <- family$linkinv
    linkinv.1  <- family$mu.eta           # 1st deriv of inverse link
    # 2nd derivative
    if (family$link=="logit") {
      linkinv.2 <- function(z) linkinv.1(z)*(1-2*linkinv(z))
    } else if (family$link=="probit") {
      linkinv.2 <- function(z) (-z)*linkinv.1(z)
    } else if (family$link=="identity") {
      linkinv.2 <- function(z) 0
    } else if (family$link=="log") {
      linkinv.2 <- function(z) linkinv(z)
    } else if (family$link=="inverse") {
      linkinv.2 <- function(z) 2*z^(-3)
    } else if (family$link=="1/mu^2") {
      linkinv.2 <- function(z) 0.75*z^(-2.5)
    } else  if (family$link=="sqrt") {
      linkinv.2 <- function(z) 2
    } else  {
      print("The specified link not allowed for computing polynomial-based CIs.")
      stop()
    }
  }


  # define the estimation method
  if (estmethod=="reg") {
    estmethod.name <- "least squares regression"
    if (is.null(vce)) vce <- "HC1"
    vce.select <- vce
  } else if (estmethod=="qreg") {
    estmethod.name <- "quantile regression"
    if (is.null(vce)) vce <- "nid"
    if (vce=="iid") {
      vce.select <- "const"
    } else {
      vce.select <- "HC1"
    }
    if (is.null(quantile)) quantile <- 0.5
  } else if (estmethod=="glm") {
    estmethod.name <- "generalized linear model"
    if (is.null(vce)) vce <- "HC1"
    vce.select <- vce
  } else {
    print("estmethod incorrectly specified.")
    stop()
  }

  ##################################error checking
  exit <- 0
  if (!is.character(binspos)) {
    if (min(binspos)<=xmin|max(binspos)>=xmax) {
      print("knots out of allowed range")
      exit <- 1
    }
  } else {
    if (binspos!="es"&binspos!="qs") {
      print("binspos incorrectly specified.")
      exit <- 1
    }
  }
  if (length(bins)==2) if (bins[1]<bins[2]) {
    print("p<s not allowed.")
    exit <- 1
  }
  if (length(testshape)==2) if (testshape[1]<testshape[2]) {
    print("p<s not allowed.")
    exit <- 1
  }
  if (length(testmodel)==2) if (testmodel[1]<testmodel[2]) {
    print("p<s not allowed.")
    exit <- 1
  }
  if (testshape[1]<deriv|testmodel[1]<deriv) {
    print("p<deriv not allowed.")
    exit <- 1
  }
  if (binsmethod!="dpi" & binsmethod!="rot") {
    print("bin selection method incorrectly specified.")
    exit <- 1
  }
  if (masspoints == "veryfew") {
    print("veryfew not allowed for testing.")
    exit <- 1
  }
  if ((length(testshape)>=1)&(length(bins)>=1)) if (testshape[1]<=bins[1]) {
    warning("p for testing <= p for bin selection not suggested.")
  }
  if ((length(testmodel)>=1)&(length(bins)>=1)) if (testmodel[1]<=bins[1]) {
    warning("p for testing <= p for bin selection not suggested.")
  }
  if (exit>0) stop()

  #################################################
  localcheck <- massadj <- T
  if (masspoints=="on") {
    localcheck <- T; massadj <- T
  } else if (masspoints=="off") {
    localcheck <- F; massadj <- F
  } else if (masspoints=="noadjust") {
    localcheck <- T; massadj <- F
  } else if (masspoints=="nolocalcheck") {
    localcheck <- F; massadj <- T
  }

  # effective size
  eN <- N <- length(x)
  Ndist <- NA
  if (massadj) {
    if (!is.null(numdist)) {
      Ndist <- numdist
    } else {
      Ndist <- length(unique(x))
    }
    eN <- min(eN, Ndist)
  }
  Nclust <- NA
  if (!is.null(cluster)) {
    if (!is.null(numclust)) {
      Nclust <- numclust
    } else {
      Nclust <- length(unique(cluster))
    }
    eN <- min(eN, Nclust)
  }

  # prepare params
  bins.p <- bins[1]; bins.s <- bins[2]
  tsha.p <- testshape[1]; tsha.s <- testshape[2]
  tmod.p <- testmodel[1]; tmod.s <- testmodel[2]
  nL <- length(testshapel); nR <- length(testshaper); nT <- length(testshape2)
  ntestshape <- nL + nR + nT

  if (binsmethod=="dpi") {
    selectmethod <- "IMSE direct plug-in"
  } else {
    selectmethod <- "IMSE rule-of-thumb"
  }

  knot <- NULL
  if (!is.character(binspos)) {
    nbins <- length(binspos)+1
    knot <- c(xmin, binspos, xmax)
    position <- "User-specified"
    es <- F
    selectmethod <- "User-specified"
  } else {
    if (binspos == "es") {
      es <- T
      position <- "Evenly-spaced"
    } else {
      es <- F
      position <- "Quantile-spaced"
    }
  }

  #### bin selection if needed ########
  if (is.null(nbins)) {
    # check if rot can be implemented
    if (is.null(nbinsrot)) {
      if (eN <= dfcheck[1]+bins.p+1+qrot) {
        print("too small effective sample size for bin selection.")
        stop()
      }
    }
    if (is.na(Ndist))  {
      Ndist.sel <- NULL
    } else {
      Ndist.sel <- Ndist
    }
    if (is.na(Nclust)) {
      Nclust.sel <- NULL
    } else {
      Nclust.sel <- Nclust
    }
    binselect <- binsregselect(y, x, w, deriv=deriv,
                               bins=bins, binspos=binspos,
                               binsmethod=binsmethod, nbinsrot=nbinsrot,
                               vce=vce.select, cluster=cluster, randcut=randcut,
                               dfcheck=dfcheck, masspoints=masspoints, weights=weights,
                               numdist=Ndist.sel, numclust=Nclust.sel)

    if (is.na(binselect$nbinsrot.regul)) {
      print("bin selection fails.")
      stop()
    }
    if (binsmethod == "rot") {
      nbins <- binselect$nbinsrot.regul
    } else if (binsmethod == "dpi") {
      nbins <- binselect$nbinsdpi
      if (is.na(nbins)) {
        warning("DPI selection fails. ROT choice used.")
        nbins <- binselect$nbinsrot.regul
      }
    }
  }

  # check eff size for testing
  tsha.fewobs <- F
  if (ntestshape!=0) if ((nbins-1)*(tsha.p-tsha.s+1)+tsha.p+1+dfcheck[2]>=eN) {
    tsha.fewobs <- T
    warning("too small eff. sample size for testing shape.")
  }

  tmod.fewobs <- F
  if (!is.null(testmodelparfit)|!is.null(testmodelpoly)) if ((nbins-1)*(tmod.p-tmod.s+1)+tmod.p+1+dfcheck[2]>=eN) {
    tmod.fewobs <- T
    warning("too small eff. sample size for testing models.")
  }

  # Generate knot
  if (is.null(knot)) {
    if (es) {
      knot <- genKnot.es(xmin, xmax, nbins)
    } else {
      knot <- genKnot.qs(x, nbins)
    }
  }
  # check local mass points
  knot <- c(knot[1], unique(knot[-1]))
  if (nbins!=length(knot)-1) {
    warning("repeated knots. Some bins dropped.")
    nbins <- length(knot)-1
  }
  if (localcheck) {
    uniqmin <- binsreg.checklocalmass(x, nbins, es, knot=knot) # mimic STATA
    if (ntestshape != 0) {
      if (uniqmin < tsha.p+1) {
        tsha.fewobs <- T
        warning("some bins have too few distinct values of x for testing shape.")
      }
    }
    if (!is.null(testmodelparfit)|!is.null(testmodelpoly)) {
      if (uniqmin < tmod.p+1) {
        tmod.fewobs <- T
        warning("some bins have too few distinct values of x for testing models.")
      }
    }
  }

  # adjust w variables
  if (!is.null(w)) {
    if (is.character(at)) {
      if (at=="mean") {
        eval.w <- colWeightedMeans(x=w, w=weights)
        if (!is.null(w.factor)) eval.w[w.factor] <- 0
      } else if (at=="median") {
        eval.w <- colWeightedMedians(x=w, w=weights)
        if (!is.null(w.factor)) eval.w[w.factor] <- 0
      } else if  (at=="zero") {
        eval.w <- rep(0, nwvar)
      }
    } else if (is.vector(at)) {
      eval.w <- at
    } else if (is.data.frame(at)) {
      eval.w <- eval.w
    }
  } else {
    eval.w <- NULL
  }

  # seed
  if (!is.null(simsseed)) set.seed(simsseed)
  # prepare grid for uniform inference
  x.grid <- binsreg.grid(knot, simsgrid)$eval

  #####################################
  ####### Shape restriction test ######
  #####################################
  stat.shapeL <- pval.shapeL <- NA
  stat.shapeR <- pval.shapeR <- NA
  stat.shape2 <- pval.shape2 <- NA
  if (nL>0) {
    stat.shapeL <- matrix(NA, nL, 2)
    pval.shapeL <- c()
  }
  if (nR>0) {
    stat.shapeR <- matrix(NA, nR, 2)
    pval.shapeR <- c()
  }
  if (nT>0) {
    stat.shape2 <- matrix(NA, nT, 2)
    pval.shape2 <- c()
  }

  if (ntestshape != 0 & !tsha.fewobs) {
    B    <- binsreg.spdes(eval=x, p=tsha.p, s=tsha.s, knot=knot, deriv=0)
    k    <- ncol(B)
    P    <- cbind(B, w)
    if (estmethod=="reg") {
      model <- lm(y ~ P-1, weights=weights)
    } else if (estmethod=="qreg") {
      model <- rq(y ~ P-1, tau=quantile, weights=weights)
    } else if (estmethod=="glm") {
      model <- glm(y ~ P-1, family=family, weights=weights)
    }
    beta <- model$coeff[1:k]
    basis.sha <- binsreg.spdes(eval=x.grid, p=tsha.p, s=tsha.s, knot=knot, deriv=deriv)

    if (estmethod=="glm" & (!nolink)) {
      pred.sha <- binsreg.pred(X=basis.sha, model=model, type="all",
                               vce=vce, cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
      basis.0     <- binsreg.spdes(eval=x.grid, p=tsha.p, s=tsha.s, knot=knot, deriv=0)
      fit.0       <- binsreg.pred(basis.0, model, type = "xb", vce=vce, cluster=cluster, deriv=0, wvec=eval.w)$fit
      pred.sha.0   <- linkinv.1(fit.0)

      if (asyvar | deriv==0) {
        pred.sha$se  <- pred.sha.0 * pred.sha$se
        if (deriv == 0) pred.sha$fit <- linkinv(pred.sha$fit)
        if (deriv == 1) pred.sha$fit <- pred.sha.0 * pred.sha$fit
      } else {
        basis.sha.1 <- basis.sha
        if (!is.null(eval.w)) {
          basis.sha.0 <- cbind(basis.0, outer(rep(1, nrow(basis.0)), eval.w))
          basis.sha.1 <- cbind(basis.sha.1, outer(rep(1, nrow(basis.sha.1)), rep(0, nwvar)))
        }
        basis.all <- linkinv.2(fit.0)*pred.sha$fit*basis.sha.0 + pred.sha.0*basis.sha.1
        pred.sha$fit <- pred.sha.0 * pred.sha$fit
        pred.sha$se  <- binsreg.pred(basis.all, model=model, type="se",
                                     vce=vce, cluster=cluster, avar=T)$se
      }
    } else {
      if (estmethod=="qreg") {
        pred.sha  <- binsreg.pred(basis.sha, model, type = "all", vce=vce,
                                  cluster=cluster, deriv=deriv, wvec=eval.w,
                                  is.qreg=TRUE, avar=asyvar, ...)
      } else {
        pred.sha  <- binsreg.pred(basis.sha, model, type = "all", vce=vce,
                                  cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
      }
    }
    fit.sha <- pred.sha$fit
    se.sha  <- pred.sha$se

    pos   <- !is.na(beta)
    k.new <- sum(pos)
    if (estmethod=="qreg") vcv.sha <- binsreg.vcov(model, type=vce, cluster=cluster, is.qreg=TRUE, ...)[1:k.new, 1:k.new]
    else                   vcv.sha <- binsreg.vcov(model, type=vce, cluster=cluster)[1:k.new, 1:k.new]

    for (j in 1:ntestshape) {
      if (j <= nL) {
        stat.shapeL[j,2] <- 1
        stat.shapeL[j,1] <- max((fit.sha - testshapel[j]) / se.sha)
      } else if (j <= nL+nR) {
        stat.shapeR[j-nL,2] <- 2
        stat.shapeR[j-nL,1] <- min((fit.sha - testshaper[j-nL]) / se.sha)
      } else {
        stat.shape2[j-nL-nR,2] <- 3
        if (is.infinite(lp)) {
          stat.shape2[j-nL-nR,1] <- max(abs((fit.sha - testshape2[j-nL-nR]) / se.sha))
        } else {
          stat.shape2[j-nL-nR,1] <- mean(abs((fit.sha - testshape2[j-nL-nR]) / se.sha)^lp)^(1/lp)
        }
      }
    }
    stat.shape <- NULL
    if (nL > 0) {
       stat.shape <- rbind(stat.shape, stat.shapeL)
    }
    if (nR > 0) {
       stat.shape <- rbind(stat.shape, stat.shapeR)
    }
    if (nT > 0) {
       stat.shape <- rbind(stat.shape, stat.shape2)
    }

    # Compute p-val
    Sigma.root <- lssqrtm(vcv.sha)
    num        <- basis.sha[,pos,drop=F] %*% Sigma.root
    denom.sha  <- sqrt(rowSums((basis.sha[, pos, drop=F] %*% vcv.sha) * basis.sha[, pos, drop=F]))
    pval.shape <- binsreg.pval(num, denom.sha, nsims, tstat=stat.shape, side=NULL, lp=lp)$pval
    if (nL!=0) {
       stat.shapeL <- stat.shapeL[,1]
       pval.shapeL <- pval.shape[1:nL]
    }
    if (nR!=0) {
       stat.shapeR <- stat.shapeR[,1]
       pval.shapeR <- pval.shape[(nL+1):(nL+nR)]
    }
    if (nT!=0) {
       stat.shape2 <- stat.shape2[,1]
       pval.shape2 <- pval.shape[(nL+nR+1):ntestshape]
    }
  }

  ############################
  #### Specification Test ####
  ############################
  stat.poly <- pval.poly <- NA
  stat.mod <- pval.mod <- NA
  if (!is.null(testmodelpoly)) {
    stat.poly <- matrix(NA, 1,2)
    pval.poly <- c()
  }
  if (!is.null(testmodelparfit)) {
    stat.mod <- matrix(NA,ncol(testmodelparfit)-1,2)
    pval.mod <- c()
  }

  if ((!is.null(testmodelparfit)|!is.null(testmodelpoly))&!tmod.fewobs) {
    if (tmod.p==tsha.p & tmod.s==tsha.s & exists("vcv.sha")) {
      exist.mod <- T
      #beta.mod <- beta.sha
      vcv.mod  <- vcv.sha
      fit.mod  <- fit.sha
      se.mod   <- se.sha
      basis.mod <- basis.sha
      denom.mod <- denom.sha
      pos <- pos
      k.new <- k.new
      model <- model
    } else {
      exist.mod <- F
      B    <- binsreg.spdes(eval=x, p=tmod.p, s=tmod.s, knot=knot, deriv=0)
      k    <- ncol(B)
      P    <- cbind(B, w)
      if (estmethod=="reg") {
        model  <- lm(y ~ P-1, weights=weights)
      } else if (estmethod=="qreg") {
        model  <- rq(y ~ P-1, tau=quantile, weights=weights)
      } else {
        model  <- glm(y ~ P-1, family=family, weights=weights)
      }

      beta <- model$coeff[1:k]
      pos <- !is.na(beta)
      # beta.mod <- beta[pos]
      k.new <- sum(pos)
      if (estmethod=="qreg") vcv.mod <- binsreg.vcov(model, type=vce, cluster=cluster, is.qreg = TRUE, ...)[1:k.new, 1:k.new]
      else                   vcv.mod <- binsreg.vcov(model, type=vce, cluster=cluster)[1:k.new, 1:k.new]
    }

    ######################
    # Test poly reg
    if (!is.null(testmodelpoly)) {
      if (!exist.mod) {
         basis.mod <- binsreg.spdes(eval=x.grid, p=tmod.p, s=tmod.s, knot=knot, deriv=deriv)

         if (estmethod=="glm" & (!nolink)) {
           pred.mod <- binsreg.pred(X=basis.mod, model=model, type="all",
                                    vce=vce, cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
           basis.0     <- binsreg.spdes(eval=x.grid, p=tmod.p, s=tmod.s, knot=knot, deriv=0)
           fit.0       <- binsreg.pred(basis.0, model, type = "xb", vce=vce, cluster=cluster, deriv=0, wvec=eval.w)$fit
           pred.mod.0   <- linkinv.1(fit.0)

           if (asyvar | deriv==0) {
             pred.mod$se  <- pred.mod.0 * pred.mod$se
             if (deriv == 0) pred.mod$fit <- linkinv(pred.mod$fit)
             if (deriv == 1) pred.mod$fit <- pred.mod.0 * pred.mod$fit
           } else {
             basis.mod.1 <- basis.mod
             if (!is.null(eval.w)) {
               basis.mod.0 <- cbind(basis.0, outer(rep(1, nrow(basis.0)), eval.w))
               basis.mod.1 <- cbind(basis.mod.1, outer(rep(1, nrow(basis.mod.1)), rep(0, nwvar)))
             }
             basis.all <- linkinv.2(fit.0)*pred.mod$fit*basis.mod.0 + pred.mod.0*basis.mod.1
             pred.mod$fit <- pred.mod.0 * pred.mod$fit
             pred.mod$se  <- binsreg.pred(basis.all, model=model, type="se",
                                          vce=vce, cluster=cluster, avar=T)$se
           }
         } else {
           if (estmethod=="qreg") {
             pred.mod  <- binsreg.pred(basis.mod, model, type = "all", vce=vce,
                                       cluster=cluster, deriv=deriv, wvec=eval.w,
                                       is.qreg=TRUE, avar=asyvar, ...)
           } else {
             pred.mod  <- binsreg.pred(basis.mod, model, type = "all", vce=vce,
                                       cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
           }
         }
         fit.mod <- pred.mod$fit
         se.mod  <- pred.mod$se

         denom.mod <- sqrt(rowSums((basis.mod[, pos, drop=F] %*% vcv.mod) * basis.mod[, pos, drop=F]))
      }

      # Run a poly reg
      x.p <- matrix(NA, N, testmodelpoly+1)
      for (j in 1:(testmodelpoly+1))  x.p[,j] <- x^(j-1)
      P.poly <- cbind(x.p, w)
      if (estmethod=="reg") {
        model.poly <- lm(y~P.poly-1, weights=weights)
      } else if (estmethod=="qreg") {
        model.poly <- rq(y~P.poly-1, tau=quantile, weights=weights)
      } else if (estmethod=="glm") {
        model.poly <- glm(y~P.poly-1, family=family, weights=weights)
      }
      beta.poly <- model.poly$coefficients
      poly.fit <- 0
      for (j in deriv:testmodelpoly) {
          poly.fit <- poly.fit + x.grid^(j-deriv)*beta.poly[j+1]*factorial(j)/factorial(j-deriv)
      }
      if (!is.null(eval.w) & deriv==0) poly.fit <- poly.fit + sum(eval.w * beta.poly[-(1:(testmodelpoly+1))])

      stat.poly[1,2] <- 3
      if (is.infinite(lp)) stat.poly[1,1] <- max(abs((fit.mod - poly.fit) / se.mod))
      else                 stat.poly[1,1] <- mean(abs((fit.mod - poly.fit) / se.mod)^lp)^(1/lp)

      Sigma.root   <- lssqrtm(vcv.mod)
      num          <- basis.mod[,pos,drop=F] %*% Sigma.root
      pval.poly    <- binsreg.pval(num, denom.mod, nsims, tstat=stat.poly, side=NULL, lp=lp)$pval
      stat.poly    <- stat.poly[,1]
    }

    ##################################
    # Test model in another data frame
    if (!is.null(testmodelparfit)) {
       x.grid <- testmodelparfit[,1]
       basis.mod <- binsreg.spdes(eval=x.grid, p=tmod.p, s=tmod.s, knot=knot, deriv=deriv)

       if (estmethod=="glm" & (!nolink)) {
         pred.mod <- binsreg.pred(X=basis.mod, model=model, type="all",
                                  vce=vce, cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
         basis.0     <- binsreg.spdes(eval=x.grid, p=tmod.p, s=tmod.s, knot=knot, deriv=0)
         fit.0       <- binsreg.pred(basis.0, model, type = "xb", vce=vce, cluster=cluster, deriv=0, wvec=eval.w)$fit
         pred.mod.0   <- linkinv.1(fit.0)

         if (asyvar | deriv==0) {
           pred.mod$se  <- pred.mod.0 * pred.mod$se
           if (deriv == 0) pred.mod$fit <- linkinv(pred.mod$fit)
           if (deriv == 1) pred.mod$fit <- pred.mod.0 * pred.mod$fit
         } else {
           basis.mod.1 <- basis.mod
           if (!is.null(eval.w)) {
             basis.mod.0 <- cbind(basis.0, outer(rep(1, nrow(basis.0)), eval.w))
             basis.mod.1 <- cbind(basis.mod.1, outer(rep(1, nrow(basis.mod.1)), rep(0, nwvar)))
           }
           basis.all <- linkinv.2(fit.0)*pred.mod$fit*basis.mod.0 + pred.mod.0*basis.mod.1
           pred.mod$fit <- pred.mod.0 * pred.mod$fit
           pred.mod$se  <- binsreg.pred(basis.all, model=model, type="se",
                                        vce=vce, cluster=cluster, avar=T)$se
         }
       } else {
         if (estmethod=="qreg") {
           pred.mod  <- binsreg.pred(basis.mod, model, type = "all", vce=vce,
                                     cluster=cluster, deriv=deriv, wvec=eval.w,
                                     is.qreg=TRUE, avar=asyvar, ...)
         } else {
           pred.mod  <- binsreg.pred(basis.mod, model, type = "all", vce=vce,
                                     cluster=cluster, deriv=deriv, wvec=eval.w, avar=asyvar)
         }
       }
       fit.mod <- pred.mod$fit
       se.mod  <- pred.mod$se

       denom.mod <- sqrt(rowSums((basis.mod[, pos, drop=F] %*% vcv.mod) * basis.mod[, pos, drop=F]))

       for (j in 2:ncol(testmodelparfit)) {
         stat.mod[j-1,2] <- 3
         if (is.infinite(lp)) stat.mod[j-1,1] <- max(abs((fit.mod - testmodelparfit[,j]) / se.mod))
         else                 stat.mod[j-1,1] <- mean(abs((fit.mod - testmodelparfit[,j]) / se.mod)^lp)^(1/lp)
       }
       Sigma.root <- lssqrtm(vcv.mod)
       num        <- basis.mod[,pos,drop=F] %*% Sigma.root
       pval.mod <- binsreg.pval(num, denom.mod, nsims, tstat=stat.mod, side=NULL, lp=lp)$pval
       stat.mod <- stat.mod[,1]
    }
  }

  ######################
  #######output#########
  ######################
  out <- list(testshapeL=list(testvalL=testshapel, stat.shapeL=stat.shapeL, pval.shapeL=pval.shapeL),
              testshapeR=list(testvalR=testshaper, stat.shapeR=stat.shapeR, pval.shapeR=pval.shapeR),
              testshape2=list(testval2=testshape2, stat.shape2=stat.shape2, pval.shape2=pval.shape2),
              testpoly=list(testpoly=testmodelpoly, stat.poly=stat.poly, pval.poly=pval.poly),
              testmodel=list(stat.model=stat.mod, pval.model=pval.mod),
              opt = list(bins.p=bins.p, bins.s=bins.s, deriv=deriv,
                         testshape=testshape, testmodel=testmodel,
                         binspos=position, binsmethod=selectmethod,
                         n=N, Ndist=Ndist, Nclust=Nclust,
                         nbins=nbins, knot=knot, lp=lp,
                         estmethod=estmethod.name, quantile=quantile, family=family))
  out$call   <- match.call()
  class(out) <- "CCFFbinstest"
  return(out)

}

##########################################################################
#' Internal function.
#'
#' @param x Class \code{CCFFbinstest} objects.
#'
#' @keywords internal
#' @export
#'

print.CCFFbinstest <- function(x, ...) {
  cat("Call: binstest\n\n")

  cat(paste("Sample size (n)                    =  ", x$opt$n,         "\n", sep=""))
  cat(paste("# of distinct values (Ndist)       =  ", x$opt$Ndist,     "\n", sep=""))
  cat(paste("# of clusters (Nclust)             =  ", x$opt$Nclust,    "\n", sep=""))
  cat(paste("Estimation Method (estmethod)      =  ", x$opt$estmethod, "\n", sep=""))
  cat(paste("Derivative (deriv)                 =  ", x$opt$deriv,     "\n", sep=""))
  cat(paste("Bin selection:", "\n"))
  cat(paste("  Method (binsmethod)              =  ", x$opt$binsmethod,"\n", sep=""))
  cat(paste("  Placement (binspos)              =  ", x$opt$binspos,   "\n", sep=""))
  cat(paste("  degree (p)                       =  ", x$opt$bins.p,    "\n", sep=""))
  cat(paste("  smooth (s)                       =  ", x$opt$bins.s,    "\n", sep=""))
  cat(paste("  # of bins (nbins)                =  ", x$opt$nbins,     "\n", sep=""))
  cat("\n")

}

################################################################################
#' Internal function.
#'
#' @param object Class \code{CCFFbinstest} objects.
#'
#' @keywords internal
#' @export
summary.CCFFbinstest <- function(object, ...) {
  x <- object
  args <- list(...)

  cat("Call: binstest\n\n")

  cat(paste("Sample size (n)                    =  ", x$opt$n,         "\n", sep=""))
  cat(paste("# of distinct values (Ndist)       =  ", x$opt$Ndist,     "\n", sep=""))
  cat(paste("# of clusters (Nclust)             =  ", x$opt$Nclust,    "\n", sep=""))
  cat(paste("Estimation Method (estmethod)      =  ", x$opt$estmethod, "\n", sep=""))
  cat(paste("Derivative (deriv)                 =  ", x$opt$deriv,     "\n", sep=""))
  cat(paste("Bin selection:", "\n"))
  cat(paste("  Method (binsmethod)              =  ", x$opt$binsmethod,"\n", sep=""))
  cat(paste("  Placement (binspos)              =  ", x$opt$binspos,   "\n", sep=""))
  cat(paste("  degree (p)                       =  ", x$opt$bins.p,    "\n", sep=""))
  cat(paste("  smooth (s)                       =  ", x$opt$bins.s,    "\n", sep=""))
  cat(paste("  # of bins (nbins)                =  ", x$opt$nbins,     "\n", sep=""))
  cat("\n")

  if (!is.null(x$testshapeL$testvalL)|!is.null(x$testshapeR$testvalR)|!is.null(x$testshape2$testval2)) {
    cat(paste("Shape Restriction Tests:")); cat("\n")
    cat(paste("degree (p) = ", x$opt$testshape[1], sep="")); cat(paste("   "))
    cat(paste("smooth (s) = ", x$opt$testshape[2], sep="")); cat("\n")

    if (!is.null(x$testshapeL$testvalL)) {
      cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
      cat(format("H0: sup mu <=", width = 15, justify = "left"))
      cat(format("sup T",         width = 12, justify = "right"))
      cat(format("p value",       width = 12, justify = "right"))
      cat("\n")
      cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
      cat("\n")
      for (i in 1:length(x$testshapeL$testvalL)) {
        cat(format(paste(x$testshapeL$testvalL[i]),    width = 15, justify = "centre"))
        cat(format(sprintf("%3.3f", x$testshapeL$stat.shapeL[i]), width = 12, justify = "right"))
        cat(format(sprintf("%3.3f", x$testshapeL$pval.shapeL[i]), width = 12, justify = "right"))
        cat("\n")
      }
      cat(paste(rep("-", 15 + 12 + 12), collapse=""))
      cat("\n")
      cat("\n")
    }

    if (!is.null(x$testshapeR$testvalR)) {
      cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
      cat(format("H0: inf mu >=", width = 15, justify = "left"))
      cat(format("inf T",         width = 12, justify = "right"))
      cat(format("p value",       width = 12, justify = "right"))
      cat("\n")
      cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
      cat("\n")
      for (i in 1:length(x$testshapeR$testvalR)) {
        cat(format(paste(x$testshapeR$testvalR[i]),               width = 15, justify = "centre"))
        cat(format(sprintf("%3.3f", x$testshapeR$stat.shapeR[i]), width = 12, justify = "right"))
        cat(format(sprintf("%3.3f", x$testshapeR$pval.shapeR[i]), width = 12, justify = "right"))
        cat("\n")
      }
      cat(paste(rep("-", 15 + 12 + 12), collapse=""))
      cat("\n")
      cat("\n")
    }

    if (!is.null(x$testshape2$testval2)) {
      cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
      cat(format("H0: mu =",      width = 15, justify = "left"))
      if (is.infinite(x$opt$lp)) {
        cat(format("sup |T|",       width = 12, justify = "right"))
      } else {
        cat(format(paste("L", x$opt$lp, " of |T|", sep=""),  width = 12, justify = "right"))
      }
      cat(format("p value",       width = 12, justify = "right"))
      cat("\n")
      cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
      cat("\n")
      for (i in 1:length(x$testshape2$testval2)) {
        cat(format(paste(x$testshape2$testval2[i]),               width = 15, justify = "centre"))
        cat(format(sprintf("%3.3f", x$testshape2$stat.shape2[i]), width = 12, justify = "right"))
        cat(format(sprintf("%3.3f", x$testshape2$pval.shape2[i]), width = 12, justify = "right"))
        cat("\n")
      }
      cat(paste(rep("-", 15 + 12 + 12), collapse=""))
      cat("\n")
      cat("\n")
    }
  }

  if (!is.null(x$testpoly$testpoly)|!all(is.na(x$testmodel$stat.model))) {
    cat(paste("Model Specification Tests:")); cat("\n")
    cat(paste("degree (p) = ", x$opt$testmodel[1], sep="")); cat(paste("   "))
    cat(paste("smooth (s) = ", x$opt$testmodel[2], sep="")); cat("\n")

    if (!is.null(x$testpoly$testpoly)) {
      cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
      cat(format("H0: mu =",      width = 15, justify = "left"))
      if (is.infinite(x$opt$lp)) {
        cat(format("sup |T|",     width = 12, justify = "right"))
      } else {
        cat(format(paste("L", x$opt$lp, " of |T|", sep=""),  width = 12, justify = "right"))
      }
      cat(format("p value",       width = 12, justify = "right"))
      cat("\n")
      cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
      cat("\n")
      cat(format(paste("poly. p=", x$testpoly$testpoly, sep=""),  width = 15, justify = "left"))
      cat(format(sprintf("%3.3f", x$testpoly$stat.poly), width = 12, justify = "right"))
      cat(format(sprintf("%3.3f", x$testpoly$pval.poly), width = 12, justify = "right"))
      cat("\n")
      cat(paste(rep("-", 15 + 12 + 12), collapse=""))
      cat("\n")
    }

    if (!all(is.na(x$testmodel$stat.model))) {
      cat(paste(rep("=", 15 + 12 + 12), collapse="")); cat("\n")
      cat(format("H0: mu =",      width = 15, justify = "left"))
      if (is.infinite(x$opt$lp)) {
        cat(format("sup |T|",     width = 12, justify = "right"))
      } else {
        cat(format(paste("L", x$opt$lp, " of |T|", sep=""),  width = 12, justify = "right"))
      }
      cat(format("p value",       width = 12, justify = "right"))
      cat("\n")
      cat(paste(rep("-", 15 + 12 + 12), collapse="")); cat("\n")
      cat("\n")
      for (i in 1:length(x$testmodel$stat.model)) {
        cat(format(paste("model ", i, sep=""),   width = 15, justify = "centre"))
        cat(format(sprintf("%3.3f", x$testmodel$stat.model[i]),   width = 12, justify = "right"))
        cat(format(sprintf("%3.3f", x$testmodel$pval.model[i]),   width = 12, justify = "right"))
        cat("\n")
      }
      cat(paste(rep("-", 15 + 12 + 12), collapse=""))
      cat("\n")
    }
  }

}
